Trabalho séries temporais - Modelagem Emprego Formal em MG

Obetivo

Analisar e comparar o desempenho de diferentes técnicas de modelagem de séries temporais na previsão da Série de Empregos Formais em Minas Gerais, visando atender os 4 pontos abaixo.

Implementar e aplicar os modelos de Regressão em Séries Temporais, ARIMA, Suavização Exponencial e Modelo de Espaço-Estado à série de Empregos Formais em Minas Gerais.
Calcular e comparar as medidas de erro MSE (mean square erro) para cada modelo, a fim de identificar o método que apresenta o melhor desempenho preditivo.
Analisar a adequação de cada modelo à série temporal em estudo, considerando seus pressupostos e características.
Discutir os resultados obtidos e apresentar conclusões sobre a efetividade das diferentes técnicas de modelagem na previsão da série [Insira o nome da série temporal].

Em Anexo, estará todo o codigo utilizado

Tecnica de Regressao em Série Temporal

A Regressão em Séries Temporais é uma técnica estatística que estabelece uma relação linear entre uma variável dependente (a série temporal a ser prevista) e uma ou mais variáveis independentes (fatores explicativos). Através de um modelo de regressão linear, a técnica busca identificar como as variáveis independentes influenciam o comportamento da série temporal ao longo do tempo. A Regressão em Séries Temporais é uma técnica para a análise e previsão de séries temporais, mas deve ser utilizada com cautela, considerando suas suposições e limitações.

Vantagens:

Simplicidade: Fácil interpretação e implementação.
Versatilidade: Adaptável a diversas séries temporais e tipos de dados.
Transparência: Permite identificar a relação individual entre cada variável independente e a série temporal.

Desvantagens:

Suposições: Requer que a série temporal e as variáveis independentes sigam distribuições lineares e independentes.
Sensibilidade a outliers: Influenciada por valores atípicos nos dados.
Limitações na previsão: Previsões de longo prazo podem ser menos precisas.

Modelo Final da Técnica

Após 5 tentativas anteriores, o modelo escolhido foi a sexta versao, visto que apresentou um menor AIC, o melhor r quadrado ajustado, e menor termo de erro. Este modelo é composto pelo tempo, pelo quadrado do tempo (para explicar a tendencia quadrática observada na serie), pelo fator sazonal e por um termo autoregressivo de ordem 1. Com isto, tem-se os pontos abaixo:

Modelo 6 - Termo autoregressivo na serie transformada

resid_1=rep(0,nEmp) 
for(i in 2:nEmp)
  resid_1[i]=modelo5$residuals[i-1]
modelo6= lm(lnEmp ~ t + t2 + factor(sazon) + resid_1)
summary(modelo6)

## 
## Call:
## lm(formula = lnEmp ~ t + t2 + factor(sazon) + resid_1)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.005222 -0.001443  0.000167  0.001492  0.004962 
## 
## Coefficients:
##                   Estimate Std. Error  t value Pr(>|t|)    
## (Intercept)      4.504e+00  1.537e-03 2930.595  < 2e-16 ***
## t               -2.929e-04  7.098e-05   -4.127 0.000139 ***
## t2               3.774e-05  1.046e-06   36.071  < 2e-16 ***
## factor(sazon)2   2.089e-03  1.631e-03    1.281 0.206195    
## factor(sazon)3   2.176e-03  1.631e-03    1.334 0.188363    
## factor(sazon)4   8.466e-03  1.569e-03    5.397 1.86e-06 ***
## factor(sazon)5   2.190e-02  1.568e-03   13.967  < 2e-16 ***
## factor(sazon)6   3.690e-02  1.568e-03   23.540  < 2e-16 ***
## factor(sazon)7   3.945e-02  1.568e-03   25.164  < 2e-16 ***
## factor(sazon)8   3.596e-02  1.568e-03   22.933  < 2e-16 ***
## factor(sazon)9   3.002e-02  1.645e-03   18.252  < 2e-16 ***
## factor(sazon)10  2.046e-02  1.632e-03   12.536  < 2e-16 ***
## factor(sazon)11  1.624e-02  1.631e-03    9.956 1.87e-13 ***
## factor(sazon)12  1.444e-03  1.631e-03    0.885 0.380261    
## resid_1          9.718e-01  7.136e-02   13.619  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.002579 on 50 degrees of freedom
## Multiple R-squared:  0.9975, Adjusted R-squared:  0.9968 
## F-statistic:  1434 on 14 and 50 DF,  p-value: < 2.2e-16

AIC(modelo6)

## [1] -575.4465

Teste de normalidade dos residuos

plot(modelo6)

shapiro.test(modelo6$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  modelo6$residuals
## W = 0.98736, p-value = 0.7498

Não ha evidencias para rejeitar H0 (H0: os residuos sao normalmente distribuidos)

durbinWatsonTest(modelo6)

##  lag Autocorrelation D-W Statistic p-value
##    1       0.3358639      1.279857   0.004
##  Alternative hypothesis: rho != 0

Há evidencias para rejeitar-se H0 - H0: nao ha auto-correlacao de ordem 1

acf(modelo6$residuals)

Ratifica-se o valor observado no teste de Durbin-Watson - ve-se autocorrelações significante

Previsao para o Modelo de Regressao em Series Temporais

Foi removido os ultimos 9 registros da série para fazer a previsão e comparação entre previsto e realizado.

real=lnEmp[57:65]
n6=length(lnEmp)-9 # Coloca em n o tamanho da serie menos 12 observacoes 
serie=rep(0,n6) # Cria a variavel serie
for(i in 1:n6)
  serie[i]=lnEmp[i] # Coloca na variavel serie a serie transformada por LN
serie=ts(serie,start=c(1999,04),frequency=12)
t=c(1:n6)
t2=t^2

sazon1=rep.int(4:12,1)
sazon2=rep.int(1:12,3)
sazon3=rep.int(1:11,1)
sazon6 = append(sazon1, sazon2)
sazon6 = append(sazon6, sazon3)

M5=lm(serie ~ t + t2+ factor(sazon6)) 
resid_m5=rep(0,n6) 
for(i in 2:n6)
  resid_m5[i]=M5$res[i-1] 
M6= lm(serie ~ t + t2 + factor(sazon6) + resid_m5)
summary(M6)

## 
## Call:
## lm(formula = serie ~ t + t2 + factor(sazon6) + resid_m5)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0048750 -0.0012220 -0.0001147  0.0011958  0.0052264 
## 
## Coefficients:
##                   Estimate Std. Error  t value Pr(>|t|)    
## (Intercept)      4.502e+00  1.629e-03 2763.130  < 2e-16 ***
## t                3.143e-05  8.442e-05    0.372    0.712    
## t2               3.148e-05  1.437e-06   21.909  < 2e-16 ***
## factor(sazon6)2  2.185e-03  1.734e-03    1.260    0.215    
## factor(sazon6)3  2.153e-03  1.734e-03    1.242    0.221    
## factor(sazon6)4  8.771e-03  1.655e-03    5.300 4.25e-06 ***
## factor(sazon6)5  2.151e-02  1.653e-03   13.013 3.80e-16 ***
## factor(sazon6)6  3.607e-02  1.652e-03   21.829  < 2e-16 ***
## factor(sazon6)7  3.763e-02  1.652e-03   22.779  < 2e-16 ***
## factor(sazon6)8  3.310e-02  1.652e-03   20.038  < 2e-16 ***
## factor(sazon6)9  2.695e-02  1.653e-03   16.310  < 2e-16 ***
## factor(sazon6)10 2.013e-02  1.654e-03   12.173 3.38e-15 ***
## factor(sazon6)11 1.599e-02  1.656e-03    9.656 4.06e-12 ***
## factor(sazon6)12 1.968e-03  1.734e-03    1.135    0.263    
## resid_m5         8.793e-01  7.456e-02   11.794 9.34e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.002451 on 41 degrees of freedom
## Multiple R-squared:  0.9962, Adjusted R-squared:  0.9949 
## F-statistic: 765.8 on 14 and 41 DF,  p-value: < 2.2e-16

previsao = rep(1,9) 
previsao[1] = M6$coef[1] + M6$coef[2]*(n6+1) + M6$coef[3]*(n6+1)^2 + M6$coef[14] + M6$coef[15]*M5$res[n6]
previsao[2] = M6$coef[1] + M6$coef[2]*(n6+2) + M6$coef[3]*(n6+2)^2 + M6$coef[15]*M5$res[n6]


for(i in 3:9)
  previsao[i] = M6$coef[1] + M6$coef[2]*(n6+i) + M6$coef[3]*(n6+i)^2 + M6$coef[i+1] + (M6$coef[15])^i*M6$res[n6]


previsao

## [1] 4.608653 4.610338 4.615821 4.619523 4.629942 4.646553 4.665044 4.670609
## [9] 4.670153

exp(previsao)

## [1] 100.3489 100.5181 101.0708 101.4456 102.5082 104.2252 106.1703 106.7627
## [9] 106.7141

previsao=ts(previsao,start = c(2003,12), frequency = 12)
plot(exp(previsao))

real=ts(real,start = c(2003,12), frequency = 12)
plot(exp(real))

previsao

##           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
## 2003                                                                        
## 2004 4.610338 4.615821 4.619523 4.629942 4.646553 4.665044 4.670609 4.670153
##      Sep Oct Nov      Dec
## 2003             4.608653
## 2004

real

##           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
## 2003                                                                        
## 2004 4.609162 4.615121 4.620059 4.630838 4.652054 4.673763 4.685828 4.692265
##      Sep Oct Nov      Dec
## 2003             4.605170
## 2004

Intervalos de previsao

li=rep(1,9)
ls=rep(1,9)
previsaofinal=rep(1,9)
p=length(M6$coef)

PHI=1
for(i in 1:9)
{
  li[i]=previsao[i] - qt(0.975,n6-p)*sd(M6$res)*sqrt(PHI)
  ls[i]=previsao[i] + qt(0.975,n6-p)*sd(M6$res)*sqrt(PHI)
  PHI=PHI+(M6$coef[15])^2*i
}
li

## [1] 4.604379 4.604646 4.608033 4.609372 4.617311 4.631382 4.647297 4.650265
## [9] 4.647200

for(i in 1:9)
{
  previsaofinal[i] = exp(previsao[i])
  li[i]=exp(li[i])
  ls[i]=exp(ls[i])
}
length(previsaofinal)

## [1] 9

length(real)

## [1] 9

vetor_prev = cbind(exp(real),previsaofinal, li, ls)
vetor_prev

##          exp(real) previsaofinal        li       ls
## Dec 2003     100.0      100.3489  99.92087 100.7788
## Jan 2004     100.4      100.5181  99.94755 101.0919
## Feb 2004     101.0      101.0708 100.28673 101.8610
## Mar 2004     101.5      101.4456 100.42111 102.4806
## Apr 2004     102.6      102.5082 101.22151 103.8112
## May 2004     104.8      104.2252 102.65580 105.8185
## Jun 2004     107.1      106.1703 104.30270 108.0713
## Jul 2004     108.4      106.7627 104.61275 108.9568
## Aug 2004     109.1      106.7141 104.29251 109.1918

Calculo do Erro Quadratico Medio de Previsao:

for(i in 1:9)
  SQP=(exp(real[i]) - previsaofinal[i])^2/9
SQP

## [1] 0.6325209

O valor observado acima é pequeno, o que mostra uma aderencia boa para o modelo de regressão em séries temporais.

Tecnica SARIMA

O modelo SARIMA (Seasonal Autoregressive Integrated Moving Average) combina a robustez do modelo ARIMA com a flexibilidade da Suavização Exponencial para prever séries temporais com sazonalidade. Essa técnica poderosa se destaca pela capacidade de capturar padrões sazonais complexos e adaptar-se a mudanças no comportamento da série ao longo do tempo. O modelo SARIMA é utilizado para previsão de séries temporais com sazonalidade. Sua capacidade de capturar padrões sazonais complexos e se adaptar a mudanças no comportamento da série o torna uma escolha popular para diversas aplicações. No entanto, é importante considerar as suposições do modelo e a necessidade de uma seleção cuidadosa dos parâmetros para garantir a precisão das previsõ

Vantagens

Adaptabilidade: Captura padrões sazonais complexos e se adapta a mudanças no comportamento da série ao longo do tempo.
Robustez: Combina a robustez do ARIMA com a flexibilidade da Suavização Exponencial.
Precisão: Alcança alta precisão em previsões de curto e médio prazo.

Desvantagens

Seleção de Parâmetros: Requer a seleção cuidadosa dos parâmetros do modelo ARIMA e dos termos sazonais.
Suposições: Requer condições de inversibilidade e estacionariedade dos componentes MA e AR, respectivamente.
Complexidade: A implementação do modelo pode ser mais complexa do que outras técnicas.

Modelo Final da Técnica

Foram feitos alguns testes de melhor ajuste dos parametros Auto Regressivo, Diferencial, Média Movel e Sazonal. Após as analises das saidas do modelo, observa-se que o modelo de melhor ajuste para a serie do Emprego Formal de Minas Gerais, foi o modelo abaixo, que apresentou o menor AIC.

fit2 <- sarima(lnEmp,1,1,0,1,1,1,6)

## initial  value -4.129883 
## iter   2 value -5.469381
## iter   3 value -5.566620
## iter   4 value -5.571658
## iter   5 value -5.576403
## iter   6 value -5.578365
## iter   7 value -5.578528
## iter   8 value -5.578827
## iter   9 value -5.578846
## iter  10 value -5.578851
## iter  11 value -5.578854
## iter  12 value -5.578859
## iter  12 value -5.578859
## iter  12 value -5.578859
## final  value -5.578859 
## converged
## initial  value -5.475502 
## iter   2 value -5.478572
## iter   3 value -5.482324
## iter   4 value -5.488341
## iter   5 value -5.489446
## iter   6 value -5.489628
## iter   7 value -5.489638
## iter   8 value -5.489640
## iter   8 value -5.489640
## final  value -5.489640 
## converged
## <><><><><><><><><><><><><><>
##  
## Coefficients: 
##      Estimate     SE  t.value p.value
## ar1    0.5956 0.1249   4.7683  0.0000
## sar1  -0.8863 0.0539 -16.4451  0.0000
## sma1  -0.5848 0.2012  -2.9066  0.0053
## 
## sigma^2 estimated as 1.27053e-05 on 55 degrees of freedom 
##  
## AIC = -8.003472  AICc = -7.995809  BIC = -7.861373 
##

fit2

## $fit
## 
## Call:
## arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D, Q), period = S), 
##     include.mean = !no.constant, transform.pars = trans, fixed = fixed, optim.control = list(trace = trc, 
##         REPORT = 1, reltol = tol))
## 
## Coefficients:
##          ar1     sar1     sma1
##       0.5956  -0.8863  -0.5848
## s.e.  0.1249   0.0539   0.2012
## 
## sigma^2 estimated as 1.271e-05:  log likelihood = 236.1,  aic = -464.2
## 
## $degrees_of_freedom
## [1] 55
## 
## $ttable
##      Estimate     SE  t.value p.value
## ar1    0.5956 0.1249   4.7683  0.0000
## sar1  -0.8863 0.0539 -16.4451  0.0000
## sma1  -0.5848 0.2012  -2.9066  0.0053
## 
## $ICs
##       AIC      AICc       BIC 
## -8.003472 -7.995809 -7.861373

shapiro.test(fit2$fit$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  fit2$fit$residuals
## W = 0.96492, p-value = 0.06214

Nota-se que o modelo não possui media movel para componente nao sazonal, e isto devido ao fato desta componente não apresentar significancia na analise da versão anterior do modelo. Através da análise da saída do modelo, é possivel comprovar a condição de inversibilidade, para o componente MA sazonal, e também é possivel comprovar a estacionariedade dos componentes autoregressivos, seja sazonal ou não sazonal.

Previsao para o SARIMA.

Da mesma maneira que a técnica anterior, exclui-se as 9 ultimas observações e faz-se a comparação entre o real e o previsto, para obtenção do indicador de qualidade de ajuste (MSE). Importante ressaltar que a série utilizada para o ajuste do SARIMA também está transformada por logaritmo.

serie

##           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
## 1999                            4.511958 4.526127 4.539030 4.541165 4.536891
## 2000 4.503137 4.506454 4.506454 4.514151 4.528289 4.548600 4.555980 4.554929
## 2001 4.518522 4.521789 4.522875 4.530447 4.541165 4.555980 4.558079 4.553877
## 2002 4.535820 4.541165 4.546481 4.557030 4.575741 4.591071 4.594109 4.593098
## 2003 4.574711 4.578826 4.579852 4.587006 4.601162 4.619073 4.621044 4.618086
##           Sep      Oct      Nov      Dec
## 1999 4.526127 4.519612 4.517431 4.503137
## 2000 4.549657 4.537961 4.534748 4.520701
## 2001 4.549657 4.548600 4.544358 4.533674
## 2002 4.591071 4.585987 4.584967 4.574711
## 2003 4.619073 4.619073 4.619073

serieSarima <- sarima.for(serie, 9, 1, 1, 0, 1, 1, 1, 6)

serieSarima$pred

##           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
## 2003                                                                        
## 2004 4.609670 4.614122 4.613904 4.618848 4.631834 4.647129 4.649037 4.646800
##      Sep Oct Nov      Dec
## 2003             4.608695
## 2004

exp(serieSarima$pred)

##           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
## 2003                                                                        
## 2004 100.4510 100.8992 100.8772 101.3772 102.7023 104.2851 104.4843 104.2508
##      Sep Oct Nov      Dec
## 2003             100.3531
## 2004

Novamente calcula-se o erro quadrático médio

for(i in 1:9)
  SQP=(exp(real[i]) - exp(serieSarima$pred[i]))^2/9
SQP

## [1] 2.61274

E já é possivel notar que embora seja usada uma técnica mais robusta, o ajuste deste modelo SARIMA ficou pior que o ajuste do modelo de regressão em séries temporais, observando o MSE de ambos.

Técnica Suavização Exponencial

A Suavização Exponencial é uma técnica de previsão robusta e adaptável que se destaca pela simplicidade e eficiência. Através da atribuição de pesos exponenciais decrescentes às observações passadas, a técnica captura as tendências recentes da série temporal e gera previsões precisas, especialmente para séries com ruído e mudanças bruscas

Vantagens

Simplicidade: Fácil implementação e interpretação.
Adaptabilidade: Rápida adaptação a mudanças recentes na série temporal.
Robustez: Eficaz para séries com ruído e mudanças bruscas.
Eficiência Computacional: Baixo custo computacional, ideal para grandes conjuntos de dados.

Desvantagens

Suposições: Requer que a série temporal siga distribuições lineares e independentes.

Modelo Final da Técnica

De fato a implementação computacional em código R para esta técnica é muito simples, bastando apenas o uso de uma função.

fitHT <- HoltWinters(lnEmp)
plot(fitted(fitHT))

plot(fitHT)

Observando o gráfico gerado pela saída do modelo, é possivel observar os 3 componentes, um nivel, uma tendencia e uma sazonalidade, além de confirmar o bom ajuste, linha vermelha, seguindo a linha preta.

Previsão para Suaviação Exponencial

Seguindo o padrão ja realizado anteriormente neste relatório, utiliza-se o vetor com as ultimas 9 observações removidas, ajusta-se um novo modelo para fazer a comparação com o realizado.

fitHT2 <- HoltWinters(serie)
plot(fitted(fitHT2))

plot(fitHT2)

predHT2 = predict(fitHT2,9, prediction.interval = TRUE)
exp(predHT2)

##                fit      upr       lwr
## Dec 2003  99.83901 100.6652  99.01956
## Jan 2004  99.77359 100.7571  98.79969
## Feb 2004 100.18276 101.3186  99.05965
## Mar 2004 100.38318 101.6638  99.11871
## Apr 2004 101.26748 102.6992  99.85568
## May 2004 102.81448 104.4077 101.24562
## Jun 2004 104.60144 106.3628 102.86927
## Jul 2004 105.04856 106.9577 103.17351
## Aug 2004 105.10032 107.1504 103.08943

exp(real)

##        Jan   Feb   Mar   Apr   May   Jun   Jul   Aug Sep Oct Nov   Dec
## 2003                                                             100.0
## 2004 100.4 101.0 101.5 102.6 104.8 107.1 108.4 109.1

plot(exp(predHT2))

plot(exp(real))

Com o modelo ajustado, é possivel avaliar o erro médio quadratico.

for(i in 1:9)
  SQP=(exp(real[i]) - exp(predHT2[i,1]))^2/9
SQP

##     fit 
## 1.77749

E o valor observado foi menor que o da técnica SARIMA, porém ainda é maior que o erro médio quadratico encontrado usando a técnica de regressão em séries temporais.

Técnica de Espaço-Estado

O modelo de Espaço-Estado se destaca por sua flexibilidade e capacidade de representar sistemas dinâmicos complexos com alta precisão. Através de um sistema de equações diferenciais e matrizes, a técnica modela as relações entre variáveis observáveis e não observáveis do sistema, capturando a dinâmica interna e as interações entre os componentes. Essa característica o torna ideal para prever o comportamento futuro de sistemas que apresentam mudanças não lineares, sazonalidade e ruído. é importante considerar a complexidade da implementação, a necessidade de expertise e a disponibilidade de dados de alta qualidade para garantir a efetividade da técnica. É importante considerar a complexidade da implementação, a necessidade de expertise e a disponibilidade de dados de alta qualidade para garantir a efetividade da técnica.

Vantagens

Flexibilidade: Permite modelar sistemas dinâmicos complexos com relações não lineares e sazonalidade.
Precisão: Alcança alta precisão em previsões, especialmente para sistemas com dinâmicas complexas.
Incorporação de Informação Externa: Possibilita a integração de informações externas ao modelo, como dados de outras fontes ou conhecimento especializado.

Desvantagens

Complexidade: A implementação do modelo pode ser mais complexa do que outras técnicas.
Necessidade de Expertise: Requer conhecimento em modelagem de sistemas dinâmicos e métodos de estimação de parâmetros.
Dados: Exige uma quantidade significativa de dados de alta qualidade para a estimação precisa dos parâmetros.

Modelo Final da Técnica

Como colocado na sessão anterior, é complexo modelar usando sistemas de espaço-estado, e o código em R também é mais complexo. É necessário ter maior dominio e maior desenvolvimento de código. E para esta técnica, a série transformada gera erro na saida do modelo, ou seja, foi necessário utilizar a série original.

model<-function(u){
  mod<-dlmModSeas(frequency=12,dV=0,dW=c(exp(u[4]),rep(0,10)))+
    dlmModPoly(2,dV=exp(u[1]),dW=(exp(u[2:3])))
}
outmle=dlmMLE(tsEmp,parm=rep(0,4),model)
exp(outmle$par)

## [1] 1.667719e-06 3.801566e-02 1.067618e-02 3.517537e-08

mod=model(outmle$par)
outmodFil=dlmFilter(tsEmp,mod)

outF<-dlmFilter(tsEmp,mod)    #Filtering
outS<-dlmSmooth(tsEmp,mod)    #smoothing

par(mfrow=c(4,1))
ts.plot((outF$m[10:nEmp,1]))
title("Sazonality: Filtering")
ts.plot((outF$m[10:nEmp,12]))
title("Slope: Filtering")
ts.plot((outF$m[10:nEmp,13]))
title("Level: Filtering")
ts.plot((outF$f[10:nEmp]),ylab="Yt^",xlab="t")
title("Preditos")

par(mfrow=c(1,1))
qqnorm(residuals(outF,sd=FALSE))
qqline(residuals(outF,sd=FALSE))

tsdiag(outF)

shapiro.test(residuals(outF,sd=FALSE))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(outF, sd = FALSE)
## W = 0.94949, p-value = 0.009997

Tendo em vista os 4 graficos gerados pelo codigo acima, nota-se que os primeiros pontos para todos eles, estão, de certa forma, um pouco fora da escala, enquanto que os demais, vão, de certa forma mais alinhados. Uma caracteristica de adequação de sistemas dinâmicos na engenharia. Importante registar que o teste de Sharipo-wilk feito em cima dos residuos do modelo de Espaço-Estado, tem um p-valor que traz evidências para rejeitar H0, ou seja, os residuos não são normalmente distribuidos.

Previsão para Modelo de Espaço-Estado

Mesmo com a ressalva acima, dos residuos do modelo não ser normalmente distribuido, foi definido fazer a previsão usando o modelo ajustado para a série sem as ultimas 9 observações.

seriePura = exp(serie)
seriePura

##        Jan   Feb   Mar   Apr   May   Jun   Jul   Aug   Sep   Oct   Nov   Dec
## 1999                    91.1  92.4  93.6  93.8  93.4  92.4  91.8  91.6  90.3
## 2000  90.3  90.6  90.6  91.3  92.6  94.5  95.2  95.1  94.6  93.5  93.2  91.9
## 2001  91.7  92.0  92.1  92.8  93.8  95.2  95.4  95.0  94.6  94.5  94.1  93.1
## 2002  93.3  93.8  94.3  95.3  97.1  98.6  98.9  98.8  98.6  98.1  98.0  97.0
## 2003  97.0  97.4  97.5  98.2  99.6 101.4 101.6 101.3 101.4 101.4 101.4

outmle2=dlmMLE(seriePura,parm=rep(0,4),model)
exp(outmle2$par)

## [1] 1.948774e-06 4.435990e-02 3.268541e-03 1.822108e-07

mod2=model(outmle2$par)
outmodFil2=dlmFilter(seriePura,mod2)

outF2<-dlmFilter(seriePura,mod2)    #Filtering
outS2<-dlmSmooth(seriePura,mod2)    #smoothing

prev2=dlmForecast(outmodFil2,n=9)
prev2$f

##           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
## 2003                                                                        
## 2004 100.7009 101.2891 101.6691 102.6409 104.2304 106.0133 106.5496 106.4995
##      Sep Oct Nov      Dec
## 2003             100.4796
## 2004

previsaoMEB=as.numeric(prev2$f)

Execução do MSE, para comparação com os demais modelos gerados anteriormente.

for(i in 1:9)
  SQP=(exp(real[i]) - previsaoMEB[i])^2/9
SQP

## [1] 0.7514231

O valor do MSE observado pelo modelo ajustado usando a técnica de Espaço-Estado foi de 0,75, o segundo melhor gerado neste estudo.

Conclusão

Para a série temporal de empregos formais no estado de Minas Gerais entre abril de 1999 e agosto de 2004, conclui-se que o modelo usando a técnica de regressão em séries temporais, se mostrou o mais precisa. Embora o esforço para se chegar ao modelo adequado, através desta técnica, tenha sido o maior, como pode ser visto no Anexo. Em seguida, um modelo de complexa implementação de Espaço-Estado, teve o segundo melhor desepenho em precisão, seguido depois pelo modelo de suavização exponencial, o mais simples de implementar técnicamente, e por fim, o modelo SARIMA, teve o pior desempenho para prever as ultimas 9 observações com precisão. Porém, em todos os casos os desvios foram pequenos, desta forma, para está série, pensando em esforço e retorno, a técnica de suavização exponencial seria a mais adequada a ser implementada, dada a simplicidade, sem perder precisão de forma significativa.

ANEXO - Codigo completo e Analise Exploratória

Importanto bibliotecas

library(FNN)
library(corrplot)
library(ggplot2)
require(car)
library(tidyverse)
library(ks)
library(tseries)
library(moments)
library(stats)
library(forecast)
library(astsa)
library(dlm)

Definicao do ambiente de trabalho e leitura de dados

setwd("/Users/victortelles/Documents/Coursera/Especializacao - UFMG/06 - Analise de Series Temporais/Trabalho Pratico Final/serie2")

stEmp <- read.csv("SerieIEF.csv", sep=";", header = T)
nEmp = length(stEmp$yt)
ts.plot(stEmp$yt, xlab = "tempo", ylab = "Indice de emprego formal em MG")

definindo a Serie Tempotal - usando o pacote tseries

tsEmp <- ts(stEmp$yt, start = c(1999,4), frequency = 12)
tsEmp

##        Jan   Feb   Mar   Apr   May   Jun   Jul   Aug   Sep   Oct   Nov   Dec
## 1999                    91.1  92.4  93.6  93.8  93.4  92.4  91.8  91.6  90.3
## 2000  90.3  90.6  90.6  91.3  92.6  94.5  95.2  95.1  94.6  93.5  93.2  91.9
## 2001  91.7  92.0  92.1  92.8  93.8  95.2  95.4  95.0  94.6  94.5  94.1  93.1
## 2002  93.3  93.8  94.3  95.3  97.1  98.6  98.9  98.8  98.6  98.1  98.0  97.0
## 2003  97.0  97.4  97.5  98.2  99.6 101.4 101.6 101.3 101.4 101.4 101.4 100.0
## 2004 100.4 101.0 101.5 102.6 104.8 107.1 108.4 109.1

plot(tsEmp)

## Analise Exploratoria da Série Temporal

Coeficiente de Variação da Serie

cv = sd(tsEmp)/mean(tsEmp)
cv

## [1] 0.04655599

Em media a série varia em 4,7%

skewness(tsEmp)

## [1] 0.8733008

kurtosis(tsEmp)

## [1] 3.275672

Observa-se simetria em torno da média. Pela curtose, temos possibilidade de gerar bastante valores nas caldas para valores de Curtoses alto - valor de referencia 3

summary(tsEmp)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   90.30   93.10   95.10   96.35   98.90  109.10

boxplot(tsEmp)

fatorSazonal=rep(4:12,1)
fatorSazonal2=rep(1:12,4)
fatorSazonal3=rep(1:8,1)
fatorSazo = append(fatorSazonal, fatorSazonal2)
fatorSazo = append(fatorSazo, fatorSazonal3)

boxplot(tsEmp~fatorSazo,xlab="Meses",ylab="Ãndice de emprego formal em MG")

Observando o boxplot mensal, é possivel observar as diferenças pluviometricas ao longo do ano

hist(tsEmp)

Ratifica-se o valor de curtose, observando caldas longas

qqnorm(tsEmp); qqline(tsEmp)

Observando o qqplot, pode-se dizer que os dados não respeitam a distribuição normal

shapiro.test(tsEmp)

## 
##  Shapiro-Wilk normality test
## 
## data:  tsEmp
## W = 0.92401, p-value = 0.0006615

Há evidencias para rejeitar H0 (H0: de que os dados sao normalmente distribuidos)

Análise de autocorrelação

acf(tsEmp)

Observa-se que a serie tem memoria longa - alem de vermos as ondas tradicionais de sazonalidade

pacf(tsEmp)

Observa-se que a ordem 1 ultrapassa as bandas limites

Teste de memória da Série

Box.test(tsEmp, lag = 1, type="Ljung-Box")

## 
##  Box-Ljung test
## 
## data:  tsEmp
## X-squared = 55.865, df = 1, p-value = 7.76e-14

Há evidencias para rejeitar H0 (H0: o conjunto de correlação de ordem 1 é igual zero)

Box.test(tsEmp, lag = 6, type="Ljung-Box")

## 
##  Box-Ljung test
## 
## data:  tsEmp
## X-squared = 201.28, df = 6, p-value < 2.2e-16

Há evidencias para rejeitar H0 (H0: o conjunto de correlação de ordem 6 é igual zero) ou seja possui memória de mais de 6 meses

Box.test(tsEmp, lag = 12, type="Ljung-Box")

## 
##  Box-Ljung test
## 
## data:  tsEmp
## X-squared = 309.24, df = 12, p-value < 2.2e-16

Há evidencias para rejeitar H0 (H0: o conjunto de correlação de ordem 12 é igual zero) ou seja possui memória de mais de 12 meses

Box.test(tsEmp, lag = 24, type="Ljung-Box")

## 
##  Box-Ljung test
## 
## data:  tsEmp
## X-squared = 344.87, df = 24, p-value < 2.2e-16

Há evidencias para rejeitar H0 (H0: o conjunto de correlação de ordem 24 é igual zero) ou seja possui memória de mais de 24 meses

Box.test(tsEmp, lag = 36, type="Ljung-Box")

## 
##  Box-Ljung test
## 
## data:  tsEmp
## X-squared = 411.54, df = 36, p-value < 2.2e-16

Há evidencias para rejeitar H0 (H0: o conjunto de correlação de ordem 36 é igual zero) ou seja possui memória de mais de 36 meses

Inicio da Modelagem

t = c(1:nEmp)

modelo1 <- lm(tsEmp ~ t)
summary(modelo1)

## 
## Call:
## lm(formula = tsEmp ~ t)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.2653 -1.4309 -0.3491  1.0658  5.9786 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 89.37481    0.51424  173.80   <2e-16 ***
## t            0.21149    0.01355   15.61   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.049 on 63 degrees of freedom
## Multiple R-squared:  0.7946, Adjusted R-squared:  0.7913 
## F-statistic: 243.7 on 1 and 63 DF,  p-value: < 2.2e-16

AIC(modelo1)

## [1] 281.6922

plot(modelo1)

shapiro.test(modelo1$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  modelo1$residuals
## W = 0.95311, p-value = 0.01517

Há evidencias para rejeitar H0 (H0: os residuos sao normalmente distribuidos)

durbinWatsonTest(modelo1)

##  lag Autocorrelation D-W Statistic p-value
##    1       0.8390991     0.1780147       0
##  Alternative hypothesis: rho != 0

Há evidencias para rejeitar-se H0 - H0: nao ha auto-correlacao de ordem 1

acf(modelo1$residuals)

Ratifica-se o valor observado no teste de Durbin-Watson - vê-se autocorrelações significantes

Modelo 2 - Ajuste para Tendencia

t2=t^2
modelo2 <- lm(tsEmp ~ t + t2)
summary(modelo2)

## 
## Call:
## lm(formula = tsEmp ~ t + t2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.42808 -1.23632 -0.04589  1.18878  3.05655 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 92.5795009  0.5825264 158.928  < 2e-16 ***
## t           -0.0755014  0.0407283  -1.854   0.0685 .  
## t2           0.0043483  0.0005981   7.270 7.31e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.518 on 62 degrees of freedom
## Multiple R-squared:  0.8891, Adjusted R-squared:  0.8856 
## F-statistic: 248.6 on 2 and 62 DF,  p-value: < 2.2e-16

AIC(modelo2)

## [1] 243.615

plot(modelo2)

shapiro.test(modelo2$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  modelo2$residuals
## W = 0.96047, p-value = 0.03629

Há evidencias para rejeitar-se H0 - H0: nao ha auto-correlacao de ordem 1

acf(modelo2$residuals)

Ratifica-se o valor observado no teste de Durbin-Watson - vê-se autocorrelações significantes

Modelo 3 - Inclusao do componente de sazonalidade

sazon1=rep.int(4:12,1)
sazon2=rep.int(1:12,4)
sazon3=rep.int(1:8,1)
sazon = append(sazon1, sazon2)
sazon = append(sazon, sazon3)
modelo3 <- lm(tsEmp ~ t + t2 + factor(sazon))
summary(modelo3)

## 
## Call:
## lm(formula = tsEmp ~ t + t2 + factor(sazon))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.21429 -0.36498 -0.04533  0.33957  1.69858 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     90.3224275  0.3473353 260.044  < 2e-16 ***
## t               -0.0352069  0.0159769  -2.204 0.032093 *  
## t2               0.0037497  0.0002349  15.966  < 2e-16 ***
## factor(sazon)2   0.1964757  0.3688266   0.533 0.596552    
## factor(sazon)3   0.2054519  0.3688916   0.557 0.580002    
## factor(sazon)4   0.8072789  0.3547153   2.276 0.027086 *  
## factor(sazon)5   2.1062529  0.3545473   5.941 2.54e-07 ***
## factor(sazon)6   3.5810608  0.3544634  10.103 9.09e-14 ***
## factor(sazon)7   3.8483693  0.3544627  10.857 7.26e-15 ***
## factor(sazon)8   3.5248450  0.3545462   9.942 1.57e-13 ***
## factor(sazon)9   2.5991027  0.3691354   7.041 4.69e-09 ***
## factor(sazon)10  1.9455762  0.3689891   5.273 2.75e-06 ***
## factor(sazon)11  1.5445503  0.3688869   4.187 0.000112 ***
## factor(sazon)12  0.1360249  0.3688260   0.369 0.713800    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5831 on 51 degrees of freedom
## Multiple R-squared:  0.9865, Adjusted R-squared:  0.9831 
## F-statistic: 287.4 on 13 and 51 DF,  p-value: < 2.2e-16

AIC(modelo3)

## [1] 128.5813

plot(modelo3)

shapiro.test(modelo3$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  modelo3$residuals
## W = 0.97801, p-value = 0.2995

Não ha evidencias para rejeitar H0 (H0: os residuos sao normalmente distribuidos)

durbinWatsonTest(modelo3)

##  lag Autocorrelation D-W Statistic p-value
##    1       0.7918974      0.249839       0
##  Alternative hypothesis: rho != 0

há evidencias para rejeitar-se H0 - H0: nao ha auto-correlacao de ordem 1

acf(modelo3$residuals)

Ratifica-se o valor observado no teste de Durbin-Watson - vê-se autocorrelações significantes

Modelo 4 - Inclusao do termo autoregressivo de ordem 1

resid_1=rep(0,nEmp) 
for(i in 2:nEmp)
  resid_1[i]=modelo3$residuals[i-1]
modelo4= lm(tsEmp ~ t + t2 + factor(sazon) + resid_1)
summary(modelo4)

## 
## Call:
## lm(formula = tsEmp ~ t + t2 + factor(sazon) + resid_1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.52095 -0.13052  0.01346  0.14802  0.43917 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     90.4072543  0.1529018 591.276  < 2e-16 ***
## t               -0.0471847  0.0070758  -6.668 1.98e-08 ***
## t2               0.0039730  0.0001044  38.042  < 2e-16 ***
## factor(sazon)2   0.1930471  0.1622457   1.190    0.240    
## factor(sazon)3   0.1981482  0.1622749   1.221    0.228    
## factor(sazon)4   0.7854122  0.1560452   5.033 6.63e-06 ***
## factor(sazon)5   2.0822974  0.1559728  13.350  < 2e-16 ***
## factor(sazon)6   3.5545698  0.1559378  22.795  < 2e-16 ***
## factor(sazon)7   3.8188962  0.1559400  24.490  < 2e-16 ***
## factor(sazon)8   3.4919433  0.1559799  22.387  < 2e-16 ***
## factor(sazon)9   2.9559550  0.1642073  18.001  < 2e-16 ***
## factor(sazon)10  1.9531826  0.1623179  12.033 2.22e-16 ***
## factor(sazon)11  1.5500678  0.1622725   9.552 7.34e-13 ***
## factor(sazon)12  0.1390069  0.1622454   0.857    0.396    
## resid_1          1.0232201  0.0700189  14.613  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2565 on 50 degrees of freedom
## Multiple R-squared:  0.9974, Adjusted R-squared:  0.9967 
## F-statistic:  1394 on 14 and 50 DF,  p-value: < 2.2e-16

AIC(modelo4)

## [1] 22.53605

plot(modelo4)

shapiro.test(modelo4$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  modelo4$residuals
## W = 0.98257, p-value = 0.4912

Ha evidencias para rejeitar H0 (H0: os residuos sao normalmente distribuidos)

durbinWatsonTest(modelo4)

##  lag Autocorrelation D-W Statistic p-value
##    1       0.3417219      1.282271   0.002
##  Alternative hypothesis: rho != 0

Não ha auto-correlacao de ordem 1

acf(modelo4$residuals)

Ratifica-se o valor observado no teste de Durbin-Watson - ve-se autocorrelacoes significantes

modelo5 - serie Transformada

lnEmp = log(tsEmp)

modelo5= lm(lnEmp ~ t + t2 + factor(sazon))
summary(modelo5)

## 
## Call:
## lm(formula = lnEmp ~ t + t2 + factor(sazon))
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.011754 -0.003275 -0.001043  0.003426  0.013831 
## 
## Coefficients:
##                   Estimate Std. Error  t value Pr(>|t|)    
## (Intercept)      4.503e+00  3.301e-03 1364.326  < 2e-16 ***
## t               -2.003e-04  1.518e-04   -1.319   0.1929    
## t2               3.601e-05  2.232e-06   16.138  < 2e-16 ***
## factor(sazon)2   2.116e-03  3.505e-03    0.604   0.5488    
## factor(sazon)3   2.232e-03  3.505e-03    0.637   0.5271    
## factor(sazon)4   8.636e-03  3.371e-03    2.562   0.0134 *  
## factor(sazon)5   2.209e-02  3.369e-03    6.555 2.74e-08 ***
## factor(sazon)6   3.711e-02  3.368e-03   11.017 4.29e-15 ***
## factor(sazon)7   3.968e-02  3.368e-03   11.779 3.63e-16 ***
## factor(sazon)8   3.621e-02  3.369e-03   10.749 1.04e-14 ***
## factor(sazon)9   2.727e-02  3.508e-03    7.773 3.30e-10 ***
## factor(sazon)10  2.040e-02  3.506e-03    5.818 3.96e-07 ***
## factor(sazon)11  1.620e-02  3.505e-03    4.621 2.63e-05 ***
## factor(sazon)12  1.421e-03  3.505e-03    0.405   0.6869    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.005541 on 51 degrees of freedom
## Multiple R-squared:  0.9883, Adjusted R-squared:  0.9853 
## F-statistic: 331.4 on 13 and 51 DF,  p-value: < 2.2e-16

AIC(modelo5)

## [1] -476.7235

plot(modelo5)

shapiro.test(modelo5$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  modelo5$residuals
## W = 0.98798, p-value = 0.7831

Ha evidencias para rejeitar H0 (H0: os residuos sao normalmente distribuidos)

durbinWatsonTest(modelo5)

##  lag Autocorrelation D-W Statistic p-value
##    1       0.8105075     0.2567286       0
##  Alternative hypothesis: rho != 0

Há evidencias para rejeitar-se H0 - H0: nao ha auto-correlacao de ordem 1

acf(modelo5$residuals)

Ratifica-se o valor observado no teste de Durbin-Watson - ve-se autocorrelacoes significantes

Modelo 6 - Termo autoregressivo na serie transformada

resid_1=rep(0,nEmp) 
for(i in 2:nEmp)
  resid_1[i]=modelo5$residuals[i-1]
modelo6= lm(lnEmp ~ t + t2 + factor(sazon) + resid_1)
summary(modelo6)

## 
## Call:
## lm(formula = lnEmp ~ t + t2 + factor(sazon) + resid_1)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.005222 -0.001443  0.000167  0.001492  0.004962 
## 
## Coefficients:
##                   Estimate Std. Error  t value Pr(>|t|)    
## (Intercept)      4.504e+00  1.537e-03 2930.595  < 2e-16 ***
## t               -2.929e-04  7.098e-05   -4.127 0.000139 ***
## t2               3.774e-05  1.046e-06   36.071  < 2e-16 ***
## factor(sazon)2   2.089e-03  1.631e-03    1.281 0.206195    
## factor(sazon)3   2.176e-03  1.631e-03    1.334 0.188363    
## factor(sazon)4   8.466e-03  1.569e-03    5.397 1.86e-06 ***
## factor(sazon)5   2.190e-02  1.568e-03   13.967  < 2e-16 ***
## factor(sazon)6   3.690e-02  1.568e-03   23.540  < 2e-16 ***
## factor(sazon)7   3.945e-02  1.568e-03   25.164  < 2e-16 ***
## factor(sazon)8   3.596e-02  1.568e-03   22.933  < 2e-16 ***
## factor(sazon)9   3.002e-02  1.645e-03   18.252  < 2e-16 ***
## factor(sazon)10  2.046e-02  1.632e-03   12.536  < 2e-16 ***
## factor(sazon)11  1.624e-02  1.631e-03    9.956 1.87e-13 ***
## factor(sazon)12  1.444e-03  1.631e-03    0.885 0.380261    
## resid_1          9.718e-01  7.136e-02   13.619  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.002579 on 50 degrees of freedom
## Multiple R-squared:  0.9975, Adjusted R-squared:  0.9968 
## F-statistic:  1434 on 14 and 50 DF,  p-value: < 2.2e-16

AIC(modelo6)

## [1] -575.4465

plot(modelo6)

shapiro.test(modelo6$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  modelo6$residuals
## W = 0.98736, p-value = 0.7498

Não ha evidencias para rejeitar H0 (H0: os residuos sao normalmente distribuidos)

durbinWatsonTest(modelo6)

##  lag Autocorrelation D-W Statistic p-value
##    1       0.3358639      1.279857   0.004
##  Alternative hypothesis: rho != 0

Há evidencias para rejeitar-se H0 - H0: nao ha auto-correlacao de ordem 1

acf(modelo6$residuals)

Ratifica-se o valor observado no teste de Durbin-Watson - ve-se autocorrelacoes significantes

Ajuste do modelo sem as 9 ultimas observacoes

real=lnEmp[57:65]

n6=length(lnEmp)-9 
serie=rep(0,n6)
for(i in 1:n6)
  serie[i]=lnEmp[i]
serie=ts(serie,start=c(1999,04),frequency=12)
t=c(1:n6)
t2=t^2

sazon1=rep.int(4:12,1)
sazon2=rep.int(1:12,3)
sazon3=rep.int(1:11,1)
sazon6 = append(sazon1, sazon2)
sazon6 = append(sazon6, sazon3)

M5=lm(serie ~ t + t2+ factor(sazon6)) 
resid_m5=rep(0,n6) 
for(i in 2:n6)
  resid_m5[i]=M5$res[i-1] 
M6= lm(serie ~ t + t2 + factor(sazon6) + resid_m5)
summary(M6)

## 
## Call:
## lm(formula = serie ~ t + t2 + factor(sazon6) + resid_m5)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0048750 -0.0012220 -0.0001147  0.0011958  0.0052264 
## 
## Coefficients:
##                   Estimate Std. Error  t value Pr(>|t|)    
## (Intercept)      4.502e+00  1.629e-03 2763.130  < 2e-16 ***
## t                3.143e-05  8.442e-05    0.372    0.712    
## t2               3.148e-05  1.437e-06   21.909  < 2e-16 ***
## factor(sazon6)2  2.185e-03  1.734e-03    1.260    0.215    
## factor(sazon6)3  2.153e-03  1.734e-03    1.242    0.221    
## factor(sazon6)4  8.771e-03  1.655e-03    5.300 4.25e-06 ***
## factor(sazon6)5  2.151e-02  1.653e-03   13.013 3.80e-16 ***
## factor(sazon6)6  3.607e-02  1.652e-03   21.829  < 2e-16 ***
## factor(sazon6)7  3.763e-02  1.652e-03   22.779  < 2e-16 ***
## factor(sazon6)8  3.310e-02  1.652e-03   20.038  < 2e-16 ***
## factor(sazon6)9  2.695e-02  1.653e-03   16.310  < 2e-16 ***
## factor(sazon6)10 2.013e-02  1.654e-03   12.173 3.38e-15 ***
## factor(sazon6)11 1.599e-02  1.656e-03    9.656 4.06e-12 ***
## factor(sazon6)12 1.968e-03  1.734e-03    1.135    0.263    
## resid_m5         8.793e-01  7.456e-02   11.794 9.34e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.002451 on 41 degrees of freedom
## Multiple R-squared:  0.9962, Adjusted R-squared:  0.9949 
## F-statistic: 765.8 on 14 and 41 DF,  p-value: < 2.2e-16

Previsao

previsao = rep(1,9) 
previsao[1] = M6$coef[1] + M6$coef[2]*(n6+1) + M6$coef[3]*(n6+1)^2 + M6$coef[14] + M6$coef[15]*M5$res[n6]
previsao[2] = M6$coef[1] + M6$coef[2]*(n6+2) + M6$coef[3]*(n6+2)^2 + M6$coef[15]*M5$res[n6]


for(i in 3:9)
  previsao[i] = M6$coef[1] + M6$coef[2]*(n6+i) + M6$coef[3]*(n6+i)^2 + M6$coef[i+1] + (M6$coef[15])^i*M6$res[n6]


previsao

## [1] 4.608653 4.610338 4.615821 4.619523 4.629942 4.646553 4.665044 4.670609
## [9] 4.670153

exp(previsao)

## [1] 100.3489 100.5181 101.0708 101.4456 102.5082 104.2252 106.1703 106.7627
## [9] 106.7141

previsao=ts(previsao,start = c(2003,12), frequency = 12)
plot(exp(previsao))

real=ts(real,start = c(2003,12), frequency = 12)
plot(exp(real))

previsao

##           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
## 2003                                                                        
## 2004 4.610338 4.615821 4.619523 4.629942 4.646553 4.665044 4.670609 4.670153
##      Sep Oct Nov      Dec
## 2003             4.608653
## 2004

real

##           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
## 2003                                                                        
## 2004 4.609162 4.615121 4.620059 4.630838 4.652054 4.673763 4.685828 4.692265
##      Sep Oct Nov      Dec
## 2003             4.605170
## 2004

Intervalos de previsao

li=rep(1,9)
ls=rep(1,9)
previsaofinal=rep(1,9)
p=length(M6$coef)

PHI=1
for(i in 1:9)
{
  li[i]=previsao[i] - qt(0.975,n6-p)*sd(M6$res)*sqrt(PHI)
  ls[i]=previsao[i] + qt(0.975,n6-p)*sd(M6$res)*sqrt(PHI)
  PHI=PHI+(M6$coef[15])^2*i
}
li

## [1] 4.604379 4.604646 4.608033 4.609372 4.617311 4.631382 4.647297 4.650265
## [9] 4.647200

for(i in 1:9)
{
  previsaofinal[i] = exp(previsao[i])
  li[i]=exp(li[i])
  ls[i]=exp(ls[i])
}
length(previsaofinal)

## [1] 9

length(real)

## [1] 9

vetor_prev = cbind(exp(real),previsaofinal, li, ls)
vetor_prev

##          exp(real) previsaofinal        li       ls
## Dec 2003     100.0      100.3489  99.92087 100.7788
## Jan 2004     100.4      100.5181  99.94755 101.0919
## Feb 2004     101.0      101.0708 100.28673 101.8610
## Mar 2004     101.5      101.4456 100.42111 102.4806
## Apr 2004     102.6      102.5082 101.22151 103.8112
## May 2004     104.8      104.2252 102.65580 105.8185
## Jun 2004     107.1      106.1703 104.30270 108.0713
## Jul 2004     108.4      106.7627 104.61275 108.9568
## Aug 2004     109.1      106.7141 104.29251 109.1918

Calculo do Erro Quadratico Medio de Previsao

for(i in 1:9)
  SQP=(exp(real[i]) - previsaofinal[i])^2/9
SQP

## [1] 0.6325209

Tecnica ARIMA(p, q, d)

pp.test(lnEmp)

## 
##  Phillips-Perron Unit Root Test
## 
## data:  lnEmp
## Dickey-Fuller Z(alpha) = -9.6712, Truncation lag parameter = 3, p-value
## = 0.5413
## alternative hypothesis: stationary

Nao tem evidencias para rejeitar a hipotese nula H0 (H0: a série tem raiz unitária, logo é estacionária)

adf.test(lnEmp)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  lnEmp
## Dickey-Fuller = -3.626, Lag order = 3, p-value = 0.03811
## alternative hypothesis: stationary

Rejeita-se a hipÃ³tese nula H0 (H0: a série não é estacionaria)

seguindo o teste Aumentado de Dickey-Fuller, A serie é estacionaria, entao necessita da componente Diferencial do ARIMA.

pacf(lnEmp)

observando o pacf, observa-se que ultrapassa-se a banda em ordem 1 - ou seja, sera necessario um componente auto-regressivo 1

Box.test(lnEmp,lag=1,type="Ljung-Box")

## 
##  Box-Ljung test
## 
## data:  lnEmp
## X-squared = 56.314, df = 1, p-value = 6.173e-14

Box.test(lnEmp,lag=6,type="Ljung-Box")

## 
##  Box-Ljung test
## 
## data:  lnEmp
## X-squared = 205.6, df = 6, p-value < 2.2e-16

Box.test(lnEmp,lag=12,type="Ljung-Box")

## 
##  Box-Ljung test
## 
## data:  lnEmp
## X-squared = 317.18, df = 12, p-value < 2.2e-16

Box.test(lnEmp,lag=24,type="Ljung-Box")

## 
##  Box-Ljung test
## 
## data:  lnEmp
## X-squared = 354.4, df = 24, p-value < 2.2e-16

Box.test(lnEmp,lag=36,type="Ljung-Box")

## 
##  Box-Ljung test
## 
## data:  lnEmp
## X-squared = 423.84, df = 36, p-value < 2.2e-16

fit1 <- sarima(lnEmp,1,1,1,1,1,1,6)

## initial  value -4.129883 
## iter   2 value -5.283240
## iter   3 value -5.544635
## iter   4 value -5.563786
## iter   5 value -5.568725
## iter   6 value -5.578846
## iter   7 value -5.589521
## iter   8 value -5.590777
## iter   9 value -5.595660
## iter  10 value -5.596667
## iter  11 value -5.596768
## iter  12 value -5.596898
## iter  13 value -5.597348
## iter  14 value -5.597612
## iter  15 value -5.597809
## iter  16 value -5.598290
## iter  17 value -5.598463
## iter  18 value -5.599008
## iter  19 value -5.599301
## iter  20 value -5.599320
## iter  21 value -5.599350
## iter  22 value -5.599390
## iter  23 value -5.599402
## iter  24 value -5.599403
## iter  25 value -5.599403
## iter  26 value -5.599404
## iter  26 value -5.599404
## iter  26 value -5.599404
## final  value -5.599404 
## converged
## initial  value -5.488519 
## iter   2 value -5.490619
## iter   3 value -5.493312
## iter   4 value -5.496938
## iter   5 value -5.497854
## iter   6 value -5.498058
## iter   7 value -5.498195
## iter   8 value -5.498292
## iter   9 value -5.498507
## iter  10 value -5.498595
## iter  11 value -5.498616
## iter  12 value -5.498618
## iter  13 value -5.498618
## iter  14 value -5.498619
## iter  15 value -5.498619
## iter  16 value -5.498620
## iter  17 value -5.498620
## iter  18 value -5.498620
## iter  18 value -5.498620
## iter  18 value -5.498620
## final  value -5.498620 
## converged
## <><><><><><><><><><><><><><>
##  
## Coefficients: 
##      Estimate     SE  t.value p.value
## ar1    0.7912 0.1807   4.3797  0.0001
## ma1   -0.2890 0.2567  -1.1258  0.2652
## sar1  -0.9040 0.0516 -17.5229  0.0000
## sma1  -0.5544 0.2138  -2.5935  0.0122
## 
## sigma^2 estimated as 1.24632e-05 on 54 degrees of freedom 
##  
## AIC = -7.986949  AICc = -7.973936  BIC = -7.809324 
##

fit1

## $fit
## 
## Call:
## arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D, Q), period = S), 
##     include.mean = !no.constant, transform.pars = trans, fixed = fixed, optim.control = list(trace = trc, 
##         REPORT = 1, reltol = tol))
## 
## Coefficients:
##          ar1      ma1     sar1     sma1
##       0.7912  -0.2890  -0.9040  -0.5544
## s.e.  0.1807   0.2567   0.0516   0.2138
## 
## sigma^2 estimated as 1.246e-05:  log likelihood = 236.62,  aic = -463.24
## 
## $degrees_of_freedom
## [1] 54
## 
## $ttable
##      Estimate     SE  t.value p.value
## ar1    0.7912 0.1807   4.3797  0.0001
## ma1   -0.2890 0.2567  -1.1258  0.2652
## sar1  -0.9040 0.0516 -17.5229  0.0000
## sma1  -0.5544 0.2138  -2.5935  0.0122
## 
## $ICs
##       AIC      AICc       BIC 
## -7.986949 -7.973936 -7.809324

ma1 nao representativo, demais componentes respeitam as condicoes e sao representativos

Mantendo a mesma condicao de sazonalidade e removendo o componente media movel 1 (ma1) do modelo - indicado pelos resultados acima

fit2 <- sarima(lnEmp,1,1,0,1,1,1,6)

## initial  value -4.129883 
## iter   2 value -5.469381
## iter   3 value -5.566620
## iter   4 value -5.571658
## iter   5 value -5.576403
## iter   6 value -5.578365
## iter   7 value -5.578528
## iter   8 value -5.578827
## iter   9 value -5.578846
## iter  10 value -5.578851
## iter  11 value -5.578854
## iter  12 value -5.578859
## iter  12 value -5.578859
## iter  12 value -5.578859
## final  value -5.578859 
## converged
## initial  value -5.475502 
## iter   2 value -5.478572
## iter   3 value -5.482324
## iter   4 value -5.488341
## iter   5 value -5.489446
## iter   6 value -5.489628
## iter   7 value -5.489638
## iter   8 value -5.489640
## iter   8 value -5.489640
## final  value -5.489640 
## converged
## <><><><><><><><><><><><><><>
##  
## Coefficients: 
##      Estimate     SE  t.value p.value
## ar1    0.5956 0.1249   4.7683  0.0000
## sar1  -0.8863 0.0539 -16.4451  0.0000
## sma1  -0.5848 0.2012  -2.9066  0.0053
## 
## sigma^2 estimated as 1.27053e-05 on 55 degrees of freedom 
##  
## AIC = -8.003472  AICc = -7.995809  BIC = -7.861373 
##

fit2

## $fit
## 
## Call:
## arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D, Q), period = S), 
##     include.mean = !no.constant, transform.pars = trans, fixed = fixed, optim.control = list(trace = trc, 
##         REPORT = 1, reltol = tol))
## 
## Coefficients:
##          ar1     sar1     sma1
##       0.5956  -0.8863  -0.5848
## s.e.  0.1249   0.0539   0.2012
## 
## sigma^2 estimated as 1.271e-05:  log likelihood = 236.1,  aic = -464.2
## 
## $degrees_of_freedom
## [1] 55
## 
## $ttable
##      Estimate     SE  t.value p.value
## ar1    0.5956 0.1249   4.7683  0.0000
## sar1  -0.8863 0.0539 -16.4451  0.0000
## sma1  -0.5848 0.2012  -2.9066  0.0053
## 
## $ICs
##       AIC      AICc       BIC 
## -8.003472 -7.995809 -7.861373

shapiro.test(fit2$fit$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  fit2$fit$residuals
## W = 0.96492, p-value = 0.06214

Não ha evidencias para rejeitar a hipotese nula de que os residuos sao normalmente distribuidos

Observa-se melhora no AIC, aumento dos graus de liberdade e menor complexidade do modelo - este deve ser usado para simulacao

Previsão utilizando melhor modelo

Previsão

serie

##           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
## 1999                            4.511958 4.526127 4.539030 4.541165 4.536891
## 2000 4.503137 4.506454 4.506454 4.514151 4.528289 4.548600 4.555980 4.554929
## 2001 4.518522 4.521789 4.522875 4.530447 4.541165 4.555980 4.558079 4.553877
## 2002 4.535820 4.541165 4.546481 4.557030 4.575741 4.591071 4.594109 4.593098
## 2003 4.574711 4.578826 4.579852 4.587006 4.601162 4.619073 4.621044 4.618086
##           Sep      Oct      Nov      Dec
## 1999 4.526127 4.519612 4.517431 4.503137
## 2000 4.549657 4.537961 4.534748 4.520701
## 2001 4.549657 4.548600 4.544358 4.533674
## 2002 4.591071 4.585987 4.584967 4.574711
## 2003 4.619073 4.619073 4.619073

serieSarima <- sarima.for(serie, 9, 1, 1, 0, 1, 1, 1, 6)

serieSarima$pred

##           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
## 2003                                                                        
## 2004 4.609670 4.614122 4.613904 4.618848 4.631834 4.647129 4.649037 4.646800
##      Sep Oct Nov      Dec
## 2003             4.608695
## 2004

exp(serieSarima$pred)

##           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
## 2003                                                                        
## 2004 100.4510 100.8992 100.8772 101.3772 102.7023 104.2851 104.4843 104.2508
##      Sep Oct Nov      Dec
## 2003             100.3531
## 2004

###Forecasting MSE

for(i in 1:9)
  SQP=(exp(real[i]) - exp(serieSarima$pred[i]))^2/9
SQP

## [1] 2.61274

Tecnica de Alisamento Exponencial

Dado tudo que foi visto ate o momento, pode-se afirmar que a serie contem uma componente de Nivel, de Tendencia e Sazonalidade, logo

fitHT <- HoltWinters(lnEmp)
plot(fitted(fitHT))

plot(fitHT)

fitHT2 <- HoltWinters(serie)
plot(fitted(fitHT2))

plot(fitHT2)

predHT2 = predict(fitHT2,9, prediction.interval = TRUE)
exp(predHT2)

##                fit      upr       lwr
## Dec 2003  99.83901 100.6652  99.01956
## Jan 2004  99.77359 100.7571  98.79969
## Feb 2004 100.18276 101.3186  99.05965
## Mar 2004 100.38318 101.6638  99.11871
## Apr 2004 101.26748 102.6992  99.85568
## May 2004 102.81448 104.4077 101.24562
## Jun 2004 104.60144 106.3628 102.86927
## Jul 2004 105.04856 106.9577 103.17351
## Aug 2004 105.10032 107.1504 103.08943

exp(real)

##        Jan   Feb   Mar   Apr   May   Jun   Jul   Aug Sep Oct Nov   Dec
## 2003                                                             100.0
## 2004 100.4 101.0 101.5 102.6 104.8 107.1 108.4 109.1

plot(exp(predHT2))

plot(exp(real))

for(i in 1:9)
  SQP=(exp(real[i]) - exp(predHT2[i,1]))^2/9
SQP

##     fit 
## 1.77749

Tecnica de Espaço-Estados

lnEmp

##           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
## 1999                            4.511958 4.526127 4.539030 4.541165 4.536891
## 2000 4.503137 4.506454 4.506454 4.514151 4.528289 4.548600 4.555980 4.554929
## 2001 4.518522 4.521789 4.522875 4.530447 4.541165 4.555980 4.558079 4.553877
## 2002 4.535820 4.541165 4.546481 4.557030 4.575741 4.591071 4.594109 4.593098
## 2003 4.574711 4.578826 4.579852 4.587006 4.601162 4.619073 4.621044 4.618086
## 2004 4.609162 4.615121 4.620059 4.630838 4.652054 4.673763 4.685828 4.692265
##           Sep      Oct      Nov      Dec
## 1999 4.526127 4.519612 4.517431 4.503137
## 2000 4.549657 4.537961 4.534748 4.520701
## 2001 4.549657 4.548600 4.544358 4.533674
## 2002 4.591071 4.585987 4.584967 4.574711
## 2003 4.619073 4.619073 4.619073 4.605170
## 2004

model<-function(u){
  mod<-dlmModSeas(frequency=12,dV=0,dW=c(exp(u[4]),rep(0,10)))+
    dlmModPoly(2,dV=exp(u[1]),dW=(exp(u[2:3])))
}
outmle=dlmMLE(tsEmp,parm=rep(0,4),model)
exp(outmle$par)

## [1] 1.667719e-06 3.801566e-02 1.067618e-02 3.517537e-08

mod=model(outmle$par)
outmodFil=dlmFilter(tsEmp,mod)

outF<-dlmFilter(tsEmp,mod)    
outS<-dlmSmooth(tsEmp,mod)    

par(mfrow=c(4,1))
ts.plot((outF$m[10:nEmp,1]))
title("Sazonality: Filtering")
ts.plot((outF$m[10:nEmp,12]))
title("Slope: Filtering")
ts.plot((outF$m[10:nEmp,13]))
title("Level: Filtering")
ts.plot((outF$f[10:nEmp]),ylab="Yt^",xlab="t")
title("Preditos")

par(mfrow=c(1,1))
myt=matrix(NA,nEmp,2)
myt[,1]=tsEmp
myt[,2]=outF$f[1:nEmp]
ts.plot(myt,ylab="Yt^",xlab="t",col=c("black","red"))

Observa-se um descasamento nos preditos em momentos iniciais, depois o modelo ajusta para as observações

ts.plot((outS$s[10:nEmp,1]))
title("Sazonality: Smoothing")

ts.plot((outS$s[10:nEmp,12]))
title("Slope: Smoothing")

ts.plot((outS$s[10:nEmp,13]))
title("Level: Smoothing")

# Análise dos resíduos

par(mfrow=c(1,1))
qqnorm(residuals(outF,sd=FALSE))
qqline(residuals(outF,sd=FALSE))

tsdiag(outF)

shapiro.test(residuals(outF,sd=FALSE))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(outF, sd = FALSE)
## W = 0.94949, p-value = 0.009997

Ha evidencias para rejeitar H0 (H0: os residuos do modelo são normalmente distribuidos)

Previsão

seriePura = exp(serie)
seriePura

##        Jan   Feb   Mar   Apr   May   Jun   Jul   Aug   Sep   Oct   Nov   Dec
## 1999                    91.1  92.4  93.6  93.8  93.4  92.4  91.8  91.6  90.3
## 2000  90.3  90.6  90.6  91.3  92.6  94.5  95.2  95.1  94.6  93.5  93.2  91.9
## 2001  91.7  92.0  92.1  92.8  93.8  95.2  95.4  95.0  94.6  94.5  94.1  93.1
## 2002  93.3  93.8  94.3  95.3  97.1  98.6  98.9  98.8  98.6  98.1  98.0  97.0
## 2003  97.0  97.4  97.5  98.2  99.6 101.4 101.6 101.3 101.4 101.4 101.4

outmle2=dlmMLE(seriePura,parm=rep(0,4),model)
exp(outmle2$par)

## [1] 1.948774e-06 4.435990e-02 3.268541e-03 1.822108e-07

mod2=model(outmle2$par)
outmodFil2=dlmFilter(seriePura,mod2)

outF2<-dlmFilter(seriePura,mod2)    
outS2<-dlmSmooth(seriePura,mod2)    

prev2=dlmForecast(outmodFil2,n=9)
prev2$f

##           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
## 2003                                                                        
## 2004 100.7009 101.2891 101.6691 102.6409 104.2304 106.0133 106.5496 106.4995
##      Sep Oct Nov      Dec
## 2003             100.4796
## 2004

previsaoMEB=as.numeric(prev2$f)
lsMEB=previsaoMEB+1.96*sqrt(as.numeric(prev2$Q[1:9]))
liMEB=previsaoMEB-1.96*sqrt(as.numeric(prev2$Q[1:9]))

vetorPrevMEB = cbind(exp(real),previsaoMEB, liMEB, lsMEB)

Forecasting MSE:

for(i in 1:9)
  SQP=(exp(real[i]) - previsaoMEB[i])^2/9
SQP

## [1] 0.7514231

A previsão obteve bons resultados com a Media do erro quadratico menor que 1

Trabalho séries temporais - Modelagem Emprego Formal em MG

Victor Telles e Leandro Barros

2024-04-26

Obetivo

Tecnica de Regressao em Série Temporal

Vantagens:

Desvantagens:

Modelo Final da Técnica

Previsao para o Modelo de Regressao em Series Temporais

Tecnica SARIMA

Vantagens

Desvantagens

Modelo Final da Técnica

Previsao para o SARIMA.

Técnica Suavização Exponencial

Vantagens

Desvantagens

Modelo Final da Técnica

Previsão para Suaviação Exponencial

Técnica de Espaço-Estado

Vantagens

Desvantagens

Modelo Final da Técnica

Previsão para Modelo de Espaço-Estado

Conclusão

ANEXO - Codigo completo e Analise Exploratória

Importanto bibliotecas

Definicao do ambiente de trabalho e leitura de dados

definindo a Serie Tempotal - usando o pacote tseries

Coeficiente de Variação da Serie

Análise de autocorrelação

Teste de memória da Série

Inicio da Modelagem

Modelo 2 - Ajuste para Tendencia

Modelo 3 - Inclusao do componente de sazonalidade

Modelo 4 - Inclusao do termo autoregressivo de ordem 1

modelo5 - serie Transformada

Modelo 6 - Termo autoregressivo na serie transformada

Ajuste do modelo sem as 9 ultimas observacoes

Previsao

Intervalos de previsao

Calculo do Erro Quadratico Medio de Previsao

Tecnica ARIMA(p, q, d)

Previsão utilizando melhor modelo

Previsão

Tecnica de Alisamento Exponencial

Tecnica de Espaço-Estados

Previsão

Forecasting MSE: